#computation of state-dependnet utilities, third life-cycle stage
function Compute_Valfunc_3(prim::Primitives, est::Estimands, res::Results, age::Int64, i_ℓp::Int64; cfact::Int64 = 0, race::Int64 = 1)
    #unpack everything
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_3, x_grid_3, h_grid, a_max_3, e_s_grid = prim #grids
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_3, n_x_3, n_ℓp, n_h, n_e_s = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, div_chars = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, α_9, τ_s, τ_p, λ_0, λ_1, λ_2, λ_3, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, α_10 = est
    @unpack γ_5, α_μ, α_x, α_amen_1, α_amen_2, α_amen_3, λ_6 = est
    @unpack emax_3 = res
    i_a, i_ap = age-44, age-43
    dist = Normal() #normal distribution

    #yank out appropriate spousal wage process
    λ_0s = λ_0s_vec[race]
    λ_1s = λ_1s_vec[race]
    λ_2s = λ_2s_vec[race]
    λ_3s = λ_3s_vec[race]
    
    #first loop over variables that affect any boolean operators so as to limit boolean evaluations
    for i_e = 1:n_e, i_x = 1:n_x_3
        e, x =  e_grid[i_e], x_grid_3[i_x]
        i_x_p = min(i_x + 1, n_x_3)


        #FIRST: MAKE SURE THAT EXPERIENCE ISN'T TOO HIGH
        if x>(age-(18 + 4*e))
            continue
        end

        for i_ℓ = 1:n_ℓ
            ℓ = ℓ_grid[i_ℓ] #standardize

            #nab location wage fixed effect
            η_h, δ, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 5], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
            d2shore, amen_index, cooling = div_chars[i_ℓp, 6], div_chars[i_ℓp, 7], div_chars[i_ℓp, 8]

            if ℓ!=1
                η_h, δ, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 5], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
                d2shore, amen_index, cooling = div_chars[i_ℓ-1, 6], div_chars[i_ℓ-1, 7], div_chars[i_ℓ-1, 8]
            end
            α_1p = α_1/coli #adjust for costs of living
            α_9p = α_9/coli
            util_bonus = α_4 * (ℓ == 1) + α_amen_1 * d2shore + α_amen_2 * amen_index + α_amen_3 * cooling #parent bonus

            for i_μ = 1:n_μ, i_m = 1:n_m, i_p = 1:n_p, i_e_s = 1:n_e_s  #loop over remainder of state space
                μ, μs, m, p, e_s = μ_grid[i_μ], μs_grid[i_m], m_grid[i_m], p_grid[i_p], e_s_grid[i_e_s] #convert to standard notation
                xh = age-(18 + 4*e_s) #experience of husband
                w = exp(λ_0 + λ_1*e + λ_2*x + λ_3*x^2 + μ + η * (1 + λ_6*e)) #agent wages
                w_s = exp(λ_0s + λ_1s*e_s + λ_2s*xh + λ_3s*xh^2 + μs + η_h) * (m>0) #spouse wsages: 0 if unmarried
                Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_4 * (m>0) #compute moving cost

                #next-period emaxes
                n_emax_0, n_emax_1 = 0, 0
                if i_a!=n_a_3 #not in terminal age
                    emax_grid = zeros((n_ℓ+1)*2)

                    #parent location option
                    emax_grid[1] = β * emax_3[i_μ, i_e, i_m, 1, 1, i_a+1, i_x, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)
                    emax_grid[2] = β * emax_3[i_μ, i_e, i_m, 2, 1, i_a+1, i_x_p, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)

                    #moving options
                    for i_ℓc = 2:(n_ℓ)
                        emax_grid[2*i_ℓc - 1] = β * emax_3[i_μ, i_e, i_m, 1, i_ℓc, i_a+1, i_x, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        emax_grid[2*i_ℓc] = β * emax_3[i_μ, i_e, i_m, 2, i_ℓc, i_a+1, i_x_p, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                    end

                    #stay8ing option
                    emax_grid[2*n_ℓ+1] = β * emax_3[i_μ, i_e, i_m, 1, i_ℓ, i_a+1, i_x, i_ℓp, i_e_s]
                    emax_grid[2*n_ℓ+2] = β * emax_3[i_μ, i_e, i_m, 2, i_ℓ, i_a+1, i_x_p, i_ℓp, i_e_s]
                    temp0, temp1 = 0, 0

                    for i_ℓc = 1:(n_ℓ+1) #location choices (parent, move, stay)
                        temp0 += exp(emax_grid[2*i_ℓc - 1])
                        temp1 += exp(emax_grid[2*i_ℓc])
                    end

                    #wrap up
                    n_emax_0 = Base.MathConstants.eulergamma + log(temp0) #add back on normalizing quantity. Works algebraically.
                    n_emax_1 = Base.MathConstants.eulergamma + log(temp1)  #done
                end
                emax_diff = n_emax_0 - n_emax_1
                cutoff, emax = 0, 0 #another preallocation

                multiplier = -1 #thing to store that keeps track of how we add in α_3 into cutoffs
                if p==1
                    multiplier = 1
                end

                #compute cutoff
                cutoff = max(-Inf, log(max((α_2 + α_10 * e + α_μ * (i_μ-1)  + α_x * x + multiplier * α_3 + α_9p * w_s + emax_diff)/α_1p, 0))) - log(w)
                res.cutoff_3[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ℓp, i_e_s] = cutoff

                #compute emax
                emax += (α_1p * w_s + α_2 + α_10 * e + α_μ * (i_μ-1) + α_x * x  + α_3 * (p==1) + α_9p * w_s + n_emax_0) * cdf(dist, cutoff/σ_ε) #not working
                emax += (α_1p * w_s + α_3 * (p==0) + n_emax_1) * (1-cdf(dist, cutoff/σ_ε)) #working, part 1
                emax += (1-cdf(dist, (cutoff - σ_ε^2)/σ_ε)) * w * exp(0.5*σ_ε^2) * α_1p #working, part 2

                #update results vector
                res.emax_3[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ℓp, i_e_s] = emax #update vector
                res.emax_3[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ℓp, i_e_s] += util_bonus #update with parent utility boost
            end
        end
    end
end

#computation of state-dependent utilities, second life-cycle stage
function Compute_Valfunc_2(prim::Primitives, est::Estimands, res::Results, age::Int64, i_ℓp::Int64; cfact::Int64 = 0, race::Int64 = 1)
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_2, x_grid_2, ac_grid, ℓp_type_grid, h_grid, a_max_2, n_τ, e_s_grid = prim #grids
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_2, n_x_2, n_ac, n_ℓp, n_h, n_e_s = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, x_grid_3, n_x_3, div_chars = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, α_9, τ_s, τ_p, λ_0, λ_1, λ_2, λ_3, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, α_10 = est
    @unpack emax_3, emax_2 = res
    @unpack γ_5, α_μ, α_x, τ_p_2, λ_4, λ_5, α_amen_1, α_amen_2, α_amen_3, λ_6 = est
    i_a, i_ap = age-40, age-39
    dist = Normal()

    #yank out appropriate spousal wage process
    λ_0s = λ_0s_vec[race]
    λ_1s = λ_1s_vec[race]
    λ_2s = λ_2s_vec[race]
    λ_3s = λ_3s_vec[race]

    #first loop over variables that affect any boolean operators so as to limit boolean evaluations
    for i_e = 1:n_e, i_x = 1:n_x_2, i_τ = 1:n_τ
        e, x =  e_grid[i_e], x_grid_2[i_x]
        i_x_p = min(i_x + 1, n_x_2)
        i_x_p_3 = i_x + 1

        #FIRST: MAKE SURE THAT EXPERIENCE ISN'T TOO HIGH
        if x>(age-(18 + 4*e))
            continue
        end

        #first loop over variables that affect any boolean operators so as to limit boolean evaluations
        for i_ℓ = 1:n_ℓ, i_ac = 1:n_ac
            ℓ, ac = ℓ_grid[i_ℓ], ac_grid[i_ac] #standardize

            #compute utility bonus associated with location states. Also keep track of what alphas we're using
            α_1p, α_2p, α_3p = α_1, α_2, α_3
            util_bonus = α_4 * (ℓ == 1)
            if i_ac>1
                α_1p, α_2p, α_3p = α_5, α_6, α_7
                util_bonus = α_8 * (ℓ == 1)
            end
            α_2p += α_10 * e #update leisure preference based on educadtion

            i_ac_p = 1 #index of next-period child age state
            if ac>=0 && ac<4 #kid aged 0 to 3
                i_ac_p = i_ac+1 #next period: kid aged one year older
            end

            #nab location wage fixed effect
            η_h, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
            d2shore, amen_index, cooling = div_chars[i_ℓp, 6], div_chars[i_ℓp, 7], div_chars[i_ℓp, 8]

            if ℓ!=1
                η_h, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
                d2shore, amen_index, cooling = div_chars[i_ℓ-1, 6], div_chars[i_ℓ-1, 7], div_chars[i_ℓ-1, 8]
            end
            α_9p = α_9/coli
            α_1p /= coli
            util_bonus += α_amen_1 * d2shore + α_amen_2 * amen_index + α_amen_3 * cooling

            #begin main loop over stage-3 state space
            for i_μ = 1:n_μ, i_m = 1:n_m, i_p = 1:n_p, i_e_s = 1:n_e_s

                #convert to standard notation
                μ, μs, m, p, e_s = μ_grid[i_μ], μs_grid[i_m], m_grid[i_m], p_grid[i_p], e_s_grid[i_e_s]
                xh = age-(18 + 4*e_s) #experience of husband

                w = exp(λ_0 + λ_1*e + λ_2*x + λ_3*x^2  + μ + η * (1 + λ_6*e) + λ_4 * (ac==0 || ac==1) + λ_5 * (ac>1)) #agent wages
                w_s = exp(λ_0s + λ_1s*e_s + λ_2s * xh + λ_3s*xh^2 + μs + η_h) * (m>0) #spuse wsages: 0 if unmarried
                Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_3 * (ac>-1) + γ_4 * (m>0)  #compute moving cost

                τ_p_star = τ_p
                if m>0
                    τ_p_star = τ_p_2
                end

                if i_τ == 2
                    τ_p_star = 0.0
                end

                #####nab CCC
                δ = div_chars[i_ℓp, 5]
                if ℓ!=1
                    δ = div_chars[i_ℓ-1, 5]
                end

                ####check for counterfactuals                
                if cfact == 2
                    δ = 0 #national subsidy
                elseif cfact == 3 && ℓ == 1 
                    δ = 0 #local subsidy
                elseif cfact == 4
                    δ/=2
                elseif cfact == 5
                    δ -= 0.25
                end

                #futz with ccc following berlinski. higher prices for married college dcuated women
                if i_m>1 && e == 1
                    δ*=1.1
                else
                    δ*=0.9
                end

                #next-period emaxes -- now can move!!!
                n_emax_0, n_emax_1 = 0, 0 #preallocation
                if i_a!=n_a_2 #not in terminal age

                    #First: hacky process to find how to normalize grids to avoid things flying off to infinity
                    emax_grid = zeros((n_ℓ+1)*2)

                    #parent location option
                    emax_grid[1] = β * emax_2[i_μ, i_e, i_m, 1, 1, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)
                    emax_grid[2] = β * emax_2[i_μ, i_e, i_m, 2, 1, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)

                    #moving options
                    for i_ℓc = 2:(n_ℓ)
                        emax_grid[2*i_ℓc - 1] = β * emax_2[i_μ, i_e, i_m, 1, i_ℓc, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        emax_grid[2*i_ℓc] = β * emax_2[i_μ, i_e, i_m, 2, i_ℓc, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                    end

                    #stay8ing option
                    emax_grid[2*n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, 1, i_ℓ, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s]
                    emax_grid[2*n_ℓ+2] = β * emax_2[i_μ, i_e, i_m, 2, i_ℓ, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s]
                    temp0, temp1 = 0, 0

                    for i_ℓc = 1:(n_ℓ+1) #location choices (parent, move, stay)
                        temp0 += exp(emax_grid[2*i_ℓc - 1])
                        temp1 += exp(emax_grid[2*i_ℓc])
                    end

                    n_emax_0 = Base.MathConstants.eulergamma + log(temp0) #add back on normalizing quantity. Works algebraically.
                    n_emax_1 = Base.MathConstants.eulergamma + log(temp1) #done
                elseif i_a == n_a_2 #in terminal age
                    emax_grid = zeros((n_ℓ+1)*2)

                    #parent location option
                    emax_grid[1] = β * emax_3[i_μ, i_e, i_m, 1, 1, 1, i_x, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)
                    emax_grid[2] = β * emax_3[i_μ, i_e, i_m, 2, 1, 1, i_x_p_3, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)

                    #moving options
                    for i_ℓc = 2:(n_ℓ)
                        emax_grid[2*i_ℓc - 1] = β * emax_3[i_μ, i_e, i_m, 1, i_ℓc, 1, i_x, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        emax_grid[2*i_ℓc] = β * emax_3[i_μ, i_e, i_m, 2, i_ℓc, 1, i_x_p_3, i_ℓp, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                    end

                    #stay8ing option
                    emax_grid[2*n_ℓ+1] = β * emax_3[i_μ, i_e, i_m, 1, i_ℓ, 1, i_x, i_ℓp, i_e_s]
                    emax_grid[2*n_ℓ+2] = β * emax_3[i_μ, i_e, i_m, 2, i_ℓ, 1, i_x_p_3, i_ℓp, i_e_s]
                    temp0, temp1 = 0, 0

                    for i_ℓc = 1:(n_ℓ+1) #location choices (parent, move, stay)
                        temp0 += exp(emax_grid[2*i_ℓc - 1])
                        temp1 += exp(emax_grid[2*i_ℓc])
                    end

                    n_emax_0 = Base.MathConstants.eulergamma + log(temp0)  #add back on normalizing quantity. Works algebraically.
                    n_emax_1 = Base.MathConstants.eulergamma + log(temp1)  #done
                end
                emax_diff = n_emax_0 - n_emax_1
                cutoff, emax = 0, 0 #another preallocation

                multiplier = -1 #thing to store that keeps track of how we add in α_3 into cutoffs
                if p==1
                    multiplier = 1
                end

                #compute cutoff
                ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star * (ℓ == 1)) * (ac>-1) #childcare costs

                if cfact == 1 
                    ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star) * (ac>-1) #grandparent transfer happens regardless of location
                end

                cutoff = max(-Inf, log(max((α_2p + multiplier * α_3p + α_μ * (i_μ-1)  + α_x * x  + α_9p * w_s + ccc * α_1p + emax_diff)/α_1p, 0))) - log(w)
                res.cutoff_2[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] = cutoff

                #compute emax
                emax += (α_1p * w_s + α_2p + α_μ * (i_μ-1) + α_x * x  + α_3p * (p==1) + α_9p * w_s + n_emax_0) * cdf(dist, cutoff/σ_ε) #not working
                emax += (α_1p * w_s + α_3p * (p==0) - ccc * α_1p + n_emax_1) * (1-cdf(dist, cutoff/σ_ε)) #working, part 1
                emax += (1-cdf(dist, (cutoff - σ_ε^2)/σ_ε)) * w * exp(0.5*σ_ε^2) * α_1p #working, part 2

                #update results vector
                res.emax_2[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] = emax #update vector
                res.emax_2[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] += util_bonus #update with parent utility boost
            end
        end
    end
end

#computation of state-dependnet utilities, first life-cycle stage. This is the hard one.
function Compute_Valfunc_1(prim::Primitives, est::Estimands, res::Results, age::Int64, i_ℓp::Int64; cfact::Int64 = 0, race::Int64 = 1)
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_1, x_grid_1, ac_grid, f_grid, ℓp_type_grid, h_grid, n_τ, e_s_grid = prim #grids
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_1, n_x_1, n_ac, n_f, n_ℓp, n_h, n_e_s = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, div_chars, x_grid_2, n_x_2, a_max_1 = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, α_9, τ_s, τ_p, λ_0, λ_1, λ_2, λ_3, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, α_10 = est
    @unpack emax_1, emax_2 = res
    @unpack γ_5, α_μ, α_x, τ_p_2, λ_4, λ_5, α_amen_1, α_amen_2, α_amen_3, λ_6 = est
    @unpack θ_1, θ_2, θ_3, θ_4, σ_θ, θ_5 = est
    dist = Normal()
    dist_fert = Normal(0, σ_θ)

    #yank out appropriate spousal wage process
    λ_0s = λ_0s_vec[race]
    λ_1s = λ_1s_vec[race]
    λ_2s = λ_2s_vec[race]
    λ_3s = λ_3s_vec[race]

    #index of current age state
    i_a = age - 21
    n_ℓc = n_ℓ + 1 #number of location choice options; add one for "stay" option

    #first loop over variables that affect any boolean operators so as to limit boolean evaluations
    for i_e = 1:n_e, i_x = 1:n_x_1, i_τ = 1:n_τ
        e, x =  e_grid[i_e], x_grid_1[i_x]
        i_x_p = min(i_x + 1, n_x_1)
        i_x_p_2 = i_x + 1

        #FIRST: MAKE SURE THAT EXPERIENCE ISN'T TOO HIGH
        if x>(age-(18 + 4*e))
            continue
        end

        #first loop over variables that affect any boolean operators so as to limit boolean evaluations
        for i_ℓ = 1:n_ℓ, i_ac = 1:n_ac
            ℓ, ac = ℓ_grid[i_ℓ], ac_grid[i_ac]#standardize

            #utility premia and parameters based on presence of child
            α_1p, α_2p, α_3p = α_1, α_2, α_3
            util_bonus = α_4 * (ℓ == 1)
            if i_ac>1
                α_1p, α_2p, α_3p = α_5, α_6, α_7
                util_bonus = α_8 * (ℓ == 1)
            end
            α_2p += α_10 * e #update leisure preference based on educadtion

            #nab location wage fixed effect
            η_h, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
            d2shore, amen_index, cooling = div_chars[i_ℓp, 6], div_chars[i_ℓp, 7], div_chars[i_ℓp, 8]

            if ℓ!=1
                η_h, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
                d2shore, amen_index, cooling = div_chars[i_ℓ-1, 6], div_chars[i_ℓ-1, 7], div_chars[i_ℓ-1, 8]
            end
            α_9p = α_9/coli
            α_1p /= coli
            util_bonus += α_amen_1 * d2shore + α_amen_2 * amen_index + α_amen_3 * cooling

            #begin loop over remainder of stage-3 state space
            for i_μ = 1:n_μ, i_m = 1:n_m, i_p = 1:n_p, i_e_s = 1:n_e_s

                #convert to standard notation
                μ, μs, m, p, e_s = μ_grid[i_μ], μs_grid[i_m], m_grid[i_m], p_grid[i_p], e_s_grid[i_e_s]
                xh = age-(18 + 4*e_s) #experience of husband
                w = exp(λ_0 + λ_1*e + λ_2*x + λ_3*x^2 + μ + η * (1 + λ_6*e) + λ_4 * (ac==0 || ac==1) + λ_5 * (ac>1)) #agent wages
                w_s = exp(λ_0s + λ_1s*e_s + λ_2s * xh + λ_3s*xh^2 + μs + η_h) * (m>0) #spuse wsages: 0 if unmarried
                Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_3 * (ac>-1) + γ_4 * (m>0) #compute moving cost

                τ_p_star = τ_p
                if m>0
                    τ_p_star = τ_p_2
                end

                if i_τ == 2
                    τ_p_star = 0.0
                end

                #####nab CCC
                δ = div_chars[i_ℓp, 5]
                if ℓ!=1
                    δ = div_chars[i_ℓ-1, 5]
                end

                 ####check for counterfactuals                
                if cfact == 2
                    δ = 0 #national subsidy
                elseif cfact == 3 && ℓ == 1 
                    δ = 0 #local subsidy
                elseif cfact == 4
                    δ/=2
                elseif cfact == 5
                    δ -= 0.25
                end

                #futz with ccc following berlinski. higher prices for married college dcuated women
                if i_m>1 && e == 1
                    δ*=1.1
                else
                    δ*=0.9
                end

                emax_fert = zeros(2) #two emaxes for two possible fertility choices
                for i_f = 1:n_f
                    f = f_grid[i_f] #zero or one
                    trans_states = [e, m, age, e_s] #things that govern probabilities of marriage, fertility

                    i_ac_p = 1 #index of next-period child age state
                    if ac>=0 && ac<4 #kid aged  0 to 3
                        i_ac_p = i_ac+1 #next period: kid aged one year older
                    end
                    if f == 1 #check fertility status
                        i_ac_p = 2 #kid aged 0 next period if woman is pregnant. This is essentially what pregnancy means
                    end
    
                    #next-period emaxes. This gets a bit hairy.
                    n_emax_0, n_emax_1 = 0, 0 #preallocation

                    if i_a == n_a_1 #in termainal period. Marriage/fertility states not changing, but can choose to move one last time

                        #First: hacky process to find how to normalize grids to avoid things flying off to infinity
                        emax_grid = zeros((n_ℓ+1)*2)

                        #parent location option
                        emax_grid[1] = β * emax_2[i_μ, i_e, i_m, 1, 1, 1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)
                        emax_grid[2] = β * emax_2[i_μ, i_e, i_m, 2, 1, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) #* (i_ℓ!=1)

                        #moving options
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[2*i_ℓc - 1] = β * emax_2[i_μ, i_e, i_m, 1, i_ℓc, 1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                            emax_grid[2*i_ℓc] = β * emax_2[i_μ, i_e, i_m, 2, i_ℓc, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        end

                        #stay8ing option
                        emax_grid[2*n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, 1, i_ℓ, 1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s]
                        emax_grid[2*n_ℓ+2] = β * emax_2[i_μ, i_e, i_m, 2, i_ℓ, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s]
                        temp0, temp1 = 0, 0

                        for i_ℓc = 1:(n_ℓ+1) #location choices (parent, move, stay)
                            temp0 += exp(emax_grid[2*i_ℓc - 1])
                            temp1 += exp(emax_grid[2*i_ℓc])
                        end

                        #wrap up
                        n_emax_0 = Base.MathConstants.eulergamma + log(temp0) #add back on normalizing quantity. Works algebraically.
                        n_emax_1 = Base.MathConstants.eulergamma + log(temp1) #done

                    elseif i_a!=n_a_1 #not in terminal period

                        #store emax values and figure out amount by which to normalize stuff
                        emax_grid = zeros(2*(n_ℓ+1))

                        #loop over potential marriage/fertility realizations
                        for i_mp = 1:n_m, i_e_s_p = 1:n_e_s
                            mp = m_grid[i_mp] #standardize
                            e_s_p = e_s_grid[i_e_s_p]
                            prob = Compute_Prob_Marr(mp, e_s_p, trans_states, prim) #compute probability of realization. Slight bottleneck here.

                            #do the hacky normalizing thing again
                            emax_grid[1] += (β * emax_1[i_μ, i_e, i_mp, 1, 1, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s_p] -  (Δ + γ_5 * div_chars[i_ℓp, 3])) * prob
                            emax_grid[2] += (β * emax_1[i_μ, i_e, i_mp, 2, 1, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓp, 3])) * prob
                            for i_ℓc = 2:(n_ℓ)
                                emax_grid[2*i_ℓc - 1] += (β * emax_1[i_μ, i_e, i_mp, 1, i_ℓc, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])) * prob
                                emax_grid[2*i_ℓc] += (β * emax_1[i_μ, i_e, i_mp, 2, i_ℓc, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])) * prob
                            end
                            emax_grid[2*n_ℓ+1] += (β * emax_1[i_μ, i_e, i_mp, 1, i_ℓ, i_a+1, i_x, i_ac_p, i_ℓp, i_τ, i_e_s_p]) * prob
                            emax_grid[2*n_ℓ+2] += (β * emax_1[i_μ, i_e, i_mp, 2, i_ℓ, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p]) * prob
                        end
                        
                        temp0, temp1 = 0, 0

                        for i_ℓc = 1:(n_ℓ+1) #location choices (parent, move, stay)
                            temp0 += exp(emax_grid[2*i_ℓc - 1])
                            temp1 += exp(emax_grid[2*i_ℓc])
                        end

                        #wrap up
                        n_emax_0 = Base.MathConstants.eulergamma + log(temp0) #done
                        n_emax_1 = Base.MathConstants.eulergamma + log(temp1) #done
                    end
                    emax_diff = n_emax_0 - n_emax_1
                    cutoff, emax = 0, 0 #another preallocation

                    multiplier = -1 #thing to store that keeps track of how we add in α_3 into cutoffs
                    if p==1
                        multiplier = 1
                    end

                    #compute cutoff
                    ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star * (ℓ == 1)) * (ac>-1) #childcare costs

                    if cfact == 1 
                        ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star) * (ac>-1) #grandparent transfer happens regardless of location
                    end

                    cutoff = max(-Inf, log(max((α_2p + multiplier * α_3p + α_μ * (i_μ-1) + α_x * x  + α_9p * w_s + ccc * α_1p + emax_diff)/α_1p, 0))) - log(w)
                    res.cutoff_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_f, i_ℓp, i_τ, i_e_s] = cutoff

                    #compute emax
                    emax += (α_1p * w_s + α_2p + α_μ * (i_μ-1) + α_x * x + α_3p * (p==1) + α_9p * w_s + n_emax_0) * cdf(dist, cutoff/σ_ε) #not working
                    emax += (α_1p * w_s + α_3p * (p==0) - ccc * α_1p + n_emax_1) * (1-cdf(dist, cutoff/σ_ε)) #working, part 1
                    emax += (1-cdf(dist, (cutoff - σ_ε^2)/σ_ε)) * w * exp(0.5*σ_ε^2) * α_1p #working, part 2
                    emax_fert[i_f] = emax  #how happy agent expects to be next period GIVEN FERTILITY CHOICE
                end

                emax_final, cutoff_fert = 0, 0
                fert_cost = θ_1 - θ_2 * (m>0) + θ_3 * age + θ_4 * age * (m>0) - θ_5 * (e)  #cost of fertility
                cutoff_fert = emax_fert[1] - emax_fert[2] + fert_cost #indifference condition for fertility decision
                mills_θ = σ_θ * (pdf(dist, cutoff_fert/σ_θ)/(1-cdf(dist, cutoff_fert/σ_θ)))

                #final emaxes
                emax_final += cdf(dist_fert, cutoff_fert) * (emax_fert[1]) #likelihood of not conceiving and associated utility
                emax_final += (1-cdf(dist_fert, cutoff_fert)) * (emax_fert[2] - fert_cost + mills_θ) #likelihood of conceiving, associated utility
                
                #update fertility cutoffs!
                res.cutoff_1_fert[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] = cutoff_fert

                #update results vector
                res.emax_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] = emax_final #update vector
                res.emax_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s] += util_bonus #update with parent utility boost
            end
        end
    end
end
